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Abstract 



In this work, we have shown the important role of the finite-size correction 
to the work function in predicting the correct positions of the centroid of 
excess charge in positively charged simple metal clusters with different 
values (2 < < 7). For this purpose, firstly we have calculated the self- 
consistent Kohn-Sham energies of neutral and singly-ionized clusters with 
sizes 2 < < 100 in the framework of local spin-density approximation 
and stabilized jellium model (S JM) as well as simple jellium model (JM) with 
rigid jellium. Secondly, we have fitted our results to the asymptotic ionization 
formulas both with and without the size correction to the work function. The 
results of fittings show that the formula containing the size correction predict 
a correct position of the centroid inside the jellium while the other predicts a 
false position outside the jellium sphere. 
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I. INTRODUCTION 



In the SJM,§ the total energy of a cluster with jellium radius R = rl^^N* is given by 

ErIu] = Ts[n] + E,,[n] + J dr {^Mn; r] + {SvUsQ{r)}[n{r) - n+(r)] (1) 

+ {eu + WR) J dr n+(r). 

The first and second terms in the right hand side of Eq. (|I]) are the non- interacting kinetic 
and the exchange-correlation energies, respectively and 

MHr)= fdr' - "fK (2) 



All equations throughout this paper are in atomic units { h = = m = qq = 1) 
unless otherwise explicitly expressed. nj^{r) is the positive jellium background density, 
n+(r) = n9{R — r) where n = 3/47rff is the bulk density. The quantity (5f)ws isS the 
average of the difference potential over the Wigner-Seitz cell and the difference potential, 
6v, is defined as the difference between the pseudopotential of a lattice of ions and the 
electrostatic potential of the jellium positive background. 6(r) takes the value of unity 
inside the jellium background and zero, outside, and wr are the Madelung energy and 
the repulsive part of the pseudo-potential, respectively. 

In the continuum approximation, the ground-state density, n/j(r), of the neutral cluster 
is given by the solution of the Euler equation! 



6n{r) 



5n(r 

n=nji ^ 



(3) 

n=nii 



In Eq (^), fi{R) is the chemical potential of electrons which is constant throughout the 
cluster for the exact ground-state density ur, and E^xc = Tg + E^c- The ground-state 
density satisfies the costraint / dr nR{r) = N*. Removing an electron from the neutral 
cluster gives rise to a new ground-state density n^(r) = nR{r) + 6n{r) with / dr 6n{r) = — 1. 

The ionization energy is defined as the difference in the ground-state energies of the two 
systems with A^* and (A^* — 1) electrons which may be represented by the expansion 
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liR) = En[n'j,] - En[nR] 



(4) 



dv 



6n{r) 



On[r) + ^J arj ar ^^^^^^^^^^^ 



6n{r)6n{r') 



n=nji 



The first term in the right hand side of Eq. @ is ~fi{R) and in local density approxi- 
mation for the kinetic and exchange-correlation energies we have 



S^Enin] 



6n(r)5n{r') 



r — r 
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^ + (5(r-r') —{nskxcin)} 



(5) 



n=?ifl(r) 



Here, ekxc{n) is the sum of the bulk kinetic and exchange-correlation energies per electron 
at the density n. Choosing a spherical distribution! for 5n{r) 



5n{r) 



(6) 



47r(i? + a)2 

and interpretting {R + a) as the centroid of the excess charge, ^(r) satisfies the following 
two equations: 

e(r) 



47r(i? + a)2 



47rr^ dr 



-1 



-irrr'^ dr = -{R + a). 



/o An{R + ay 

The simplest choice for ^(r) which stisfies the above two constraints, is 



(7) 



(8) 



-6{r-R-a) 



(9) 



Using the Eqs. (|5|)-(||) in the second term of the right hand side of the Eq. @) one 
obtains 



I{R) = -fxiR) + 



+ 0{R 



-2\ 



(10) 



2{R + a) 

Since fi{R) is a constant, independent of r, for the exact ground-state density of the 
neutral cluster, n/j(r), it is possible to evaluate it at a point deep inside the cluster where 
nf>{r) = n. The density profile rtnir) and the potential for large clusters can 

be expressed as corresponding quantities for the planar surface (i.e. R = 00) plus a size 
correction^ 



nnix)=nix) + ^ + 0{R~') (11) 

Mx)=</)(x) + ^ + 0(i?-2) (12) 

where x = r — R gives the distance from the jeUium edge. At large distances from the jeUium 
edge outside the cluster, we have /(+oo) = 0(+cxd) = /i(+oo) = 0, and for a point deep 
inside the cluster {x = — oo) we have n{—oo) = n, /(— oo) = 0. Hence Eq. (|]) reads 



/TIN ,/ N h(-oo) ,^ . d , , 

H[R) = 0(-oo) H — h {dv)ws + -Q;^{n£kxc{n)\ 



(13) 

Upon inserting Eq. ([13|) into Eq. (p^Of ) and using the definition of the work function of a 
planar surface in the sjmI 



d 

W = Acf)- {6v)ws - g;^{nekxc{n)} 



(14) 

where A0 = 0(+oo) — 0(— oo) is the electrostatic dipole barrieri, one obtains 



where c = — /;,(— oo). 

Equation ( [T5| ) is also valid in the JM for which is obtained by putting ((5w)ws = 
in Eq. (|1^) and recalculating A(f) with the JM density profile. There are two methodJii 
for calculating W and c called the "Koopmans' method" and the "change in self-consistent 
field (ASCF) method" which give the same results for the correct n{x), but the ASCF 
method has the advantage that its results are less sensitive to the exact density profile n{x). 
However, in our calculations, all the quantities W, c, and a are obtained from fitting of our 



self-consistent ionization energies of the clusters to Eq. (|15|) 



II. CALCULATIONAL SCHEME AND RESULTS 

In this work, using the JM as well as the S JM with local spin-density approximation,! we 
have solved the self-consistent Kohn-Sham equations^ and obtained the energies of neutral 



and singly- ionized A^-atomic clusters of different sizes (2 < < 100). The calculations have 
been performed for different fg values of the jellium sphere (2 < < 7) in steps 0.5 bohr. 
For each fg value, we have calculated the ionization energies 

I{R) = En{N-l)-ER{N), (16) 

for 2 < < 100 with R = N'^/^fs. Then, by fitting the I{R) values to Eq. we have 
obtained the values of W, c, and a for that value. 

Figure |l](a) shows the work function obtained for our SJM data and is labeled as SJM(l) 
because, in this paper we have called the Eq. ([I5| ) as scheme 1. To compare our results with 
the bulk calculation results , SJM(bulk), we have also plotted the data from table IV of Ref. 



T2| . As is seen, our results show a good agreement with the bulk calculations results. 

In Fig. |l](b) we have compared our JM(1) results with the bulk JM results of Ref. [12. 



Here also the good agreement is obvious. 

Figure |l](c) compares the JM(1) and the SJM(l) results. These two curves coincide at 
fs ^ 4.18 where the JM is mechanically stable. The difference M/^sjm ~ W^jm is positive for 
fs < 4.18 and changes sign for > 4.18. This behavior can easily be predicted if one resorts 
to the corresponding expressions for the work functions: For fg ^ 4.18, we have (5f )ws ~ 
and therefore, the density profiles are the same which results in A0[nsjM] = A0[njM] and 
in turn, by Eq. ( p^lf ) leads to the same workfunctions. However, since the term (5f )ws iii 
Eq. ([1^ changes sign at fg = 4.18, it gives rise to the change in the sign of M^sjm — Wjm- 

In Fig. |l](d), the values of the quantum size correction c for the JM(1) and the SJM(l) 
are compared. These two plots have two intersections. The one at ~ 4.18 is obvious if 
we consider the dependenceii of the quantity c on the density profile n{x) 

/+00 
dx{x [n{x) — n6{—x)] + xf{x)}, (17) 
-oo 

/+00 r+oo 
dxf{x) = —2 dxx[n{x) — n6{—x)] (18) 
-oo J —oo 

= (19) 



The intersection at ^ 2.5 is because of the improper behavior of the properties at high 



densities in the JM ( See for example Fig. 6(a) of Ref. |T3| for the surface energy which 
becomes negative for < 2.5). 

In Fig. 0(a), we have shown the position of the centroid of excess charge, a, for the 
singly-ionized clusters as a function of fg both for the JM(1) and the SJM(l). As is seen, 
both plots predict negative values which means that the centroid for a positively charged 
cluster lies inside the jellium [ the jellium surface is taken as the origin, see Eq. (||)]. 
This is consistent with our understandings that the excess charge in a metal resides on 
the surface. On the other hand, the JM(1) and the SJM(l) have opposite behaviors: the 
former is decreasing ( increasing in absolute value) and the latter is increasing ( decreasing 
in absolute value) and intersect each other at fg ^ 4.18. The explanation of these opposite 
behaviors is straightforward: In simple JM, if we take the thickness of the jellium shell at the 
edge as A, then the volume of this shell contains a unit charge (for singly-ionized clusters) 

A(4rf)n = 1, (20) 

which gives 

By Eq. (^Tf), for a fixed R, in JM as increases, the position of the centroid moves inside 
from the edge towards the center of the jellium sphere. However, in the S JM, there exist two 
competing effects: On the one hand, as in the simple JM, because the jellium background 
is the same, increasing fs has the tendency to increase A but, on the other hand, the SJM 
Kohn-Sham effective potential contains an extra term (5f )ws which is an increasing function 
of fg ( See Table IV of Ref. |^ ). This term is negative for fg < 4.18 and hence deepens 
the effective potential ( compared to the JM ) which reduces the spill-out of the electrons, 
so that asjM < c^jm- On the other hand, for fg > 4.18 the extra term is positive which 
shallowen the effective potential ( compared to the JM ) which in turn increases the range 
of the Kohn-sham orbitals from which the electron density is calculated. The result is that 



in the SJM, the overlap of the negative and positive charges increases (the shell thickness 
decreases) with increasing fg, so that asjM > am- 

The quantity a calculated above is, of course, the centroid position for singly-ionized 
cluster. To calculate the centroid position for doubly-ionized clusters, we fitted the Kohn- 
Sham energies of doubly-ionized clusters to the Eq. (20) of Ref. ^ which is obtained from 
Taylor expansion of the energy of a -2-ply charged cluster about z = 0. The resulted values 
for a was +0.34, -0.18, -0.29 for values 2, 4, 6, respectively. Obviousely, the result for 
fs = 2 is incorrect because of its wrong sign. The predicted value for = 6 is quite 
acceptable and consistent with the above arguments. One reason for the incorrect results 
may be in the use of Taylor expansion about z = for points 2; >> 0, as in Eq. (20) of Ref. 
^. A better approximation is obtained by using the z-th ionization energy 

= Ea{N- - z) - Eg(N' -z+l) (22) 
._,go,^._, + l) + _l_, (23) 

where we have used the Taylor expansion of the ground-state energy Ef({N* —z) at N*—z+l. 
The quantity e^^{N* — z + l) is the highest occupied Kohn-Sham orbital energyS for the 
(z — l)-ply ionized cluster with radius R = fs(N*y^^. 

It is instructive to compare the behaviors of the position of the centroid of excess charge 
with the position of the image plane of a planar surface. In Fig. ^(b) we have compared 
these quantities in the JM and the SJM. As is seen, both of them are positive and intersect 
at fg ^ 4. In the JM plot we have used the Eq. (8) of Ref. |^ and for the SJM we have used 



the data from Fig. 1 of Ref. |T4[ Here, although the argument used for the finite cluster to 
calculate A is useless but, the effect of {6v) is the same as in the finite cluster case which 
leads to similar behaviors. 

Now, we show the results of fitting our self-consistent Kohn-Sham energies to Eq. (|I5|) 
but neglecting the quantum size correction, as in Eq. (16) of Ref. ^ Putting c = in Eq. 
(P^) we call the resulting equation as scheme 2. 

In Fig. ^(a) we have compared the work functions obtained in schemes 1 and 2. As is 



seen, the scheme 2 predicts lower values than the scheme 1 otherwise the shapes are the 
same. 

Figure 0(b) shows the plots of a for JM and SJM in scheme 2. As we see, both of them 
are positive and have increasing behaviors. This means that the position of the centroid of 
the excess positive charge (for the singly-ionized clusters) lies outside the jellium edge which 
is incorrect because there exists no positive charge. Therefore, inclusion of the quantum size 
effect is vital for predicting correct values of the physical quantities. 

Finally, we introduce the scheme 3 which is obtained by putting a = in Eq. (|15|) as 



Eqs. (1) of Refs. § and [T^. In Fig. ^(a) we have compared the SJM work functions in the 
schemes 1 and 3 with the bulk value£l. It is seen that scheme 3 improves the workfunction 
and reduces the distance from the bulk values at higher densities. 

The SJM fitted values of the size correction for schemes 1 and 3 are compared in Fig. 
^(b). It is seen that in scheme 3, the values are nearly constant for the metallic densities 
c ~ —?>.leV in good agreement with the result of Ref. ^ 

III. SUMMARY AND CONCLUSIONS 

We have calculated the self-consistent Kohn-Sham energies for neutral and singly-ionized 
clusters with different sizes and repeated the calculations for different fg values in the metallic 
range. The ionization energies for each has been separately fitted to Eq. (|15D, we called 
scheme 1. The results obtained for the position of the centroid of excess charge is consistent 
with the understanding that excess charges reside at the outer surface of a metal particle. 
However, for small doubly-ionized clusters, using the Taylor expansion about z = gives 
incorrect a values for clusters of high electron-density metals. A better result would be 
obtained by using a Taylor expansion of the energy of 2;-ply ionized cluster about the point 
(^-1). 

Next, in scheme 2, we have fitted the same ionization energies to Eq. (|T5|) taking c = 0. 
The results show incorrect values for the centroid position. 
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Finally, in scheme 3, we put a = in Eq. ([T5|) and fit the ionization energies, which 
results in better values for the work function at higher densities. 

To conclude, if one wishes to calculate the position of the centroid of excess charge of 
z-ply charged cluster, one should use the Taylor expansion about point [z — 1) as in Eq. 
(P3|). On the other hand, if one wishes choose one of the versions of Eq. ( PTSD to calculate the 
ionization energies ( or electron affinities) of simple metal clusters, the best choice would be 
the Eq. ([T5|) with the ASCF values for W and c, and for sufficiently large clusters {R » a) 
the scheme 3 with the ASCF values for W and c gives the same results as scheme 2. 
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FIGURES 

FIG. 1. Work functions and quantum size corrections in electron volts from scheme 1. In (a)- 
the fitted values of the work function has been compared with those of the bulk calculation for the 
SJM. (b)- the same as (a) but for the JM. (c)- compares the fitted JM and SJM work functions, 
(d)- compares the fitted quantum size corrections in the JM and the SJM. 

FIG. 2. (a)- The scheme 1 positions of the centroid of excess charge, in atomic units, as functions 

of fg for the JM and the SJM. (b)- The positions of the image planes for the planar surfaces, in 
atomic units, from bulk calculations of the JM and the SJM. 

FIG. 3. (a)- Work functions of schemes 1 and 2 in electron volts are compared, (b)- Position 
of the centroid of excess charge, in atomic units, in scheme 2. Both the SJM and the JM show 
increasing behaviors. 

FIG. 4. (a)- Work functions in electron volts for schemes 1 and 3 are compared with the bulk, 
(b)- quantum size corrections in electrov volts for the schemes 1 and 3. The scheme 3 predicts 
more or less a constant value. 
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